clc; clear;
% 对比 k=0.04 与 k=0.06 两种增益设置的影响，进行数值积分。

params.a = 0.7;
params.b = 0.8;
params.alpha = 0.25;
params.beta = 0.01;
params.delta = 0.1;
% 设置周期性输入参
params.A = 0.35;   % 振幅
params.f = 0.72;   % 频率
params.epsilon = 0.001;

% 初始条件 [x, y, z, c]
x0 = 0.2;
y0 = 0;
z0 = 0.01;
c0 = 0.1;
y_init = [x0; y0; z0; c0];

tspan = [0, 1000];   % 仿真时间区间
h =0.05;           % 步长

% 对两种不同 k 值进行仿真

% k = 0.04 情况
params1 = params;
params1.k = 0.04;
fun1 = @(t, y) hybridNeuron_AdpC(t, y, params1);
[t1, Y1] = rk4(fun1, tspan, y_init, h);

% k = 0.06 情况
params2 = params;
params2.k = 0.06;
fun2 = @(t, y) hybridNeuron_AdpC(t, y, params2);
[t2, Y2] = rk4(fun2, tspan, y_init, h);

% 能量定义： H = 0.5*x^2 + 0.5*(1/c)*y^2 + 0.5*b*z*sin(beta*z)*y
H1 = 0.5 * Y1(:,1).^2 + 0.5 * (Y1(:,2).^2)./Y1(:,4) + 0.5 * params1.b * Y1(:,3) .* sin(params1.beta*Y1(:,3)) .* Y1(:,2);
H2 = 0.5 * Y2(:,1).^2 + 0.5 * (Y2(:,2).^2)./Y2(:,4) + 0.5 * params2.b * Y2(:,3) .* sin(params2.beta*Y2(:,3)) .* Y2(:,2);

subplot(2, 1, 1);
plot(t1, H1, 'b');
xlabel('时间 \tau');
ylabel('能量 H');

subplot(2, 1, 2);
plot(t2, H2, 'r'); 
xlabel('时间 \tau');
ylabel('能量 H');